# Monte Carlo Experiment I: EPBR Data
rm(list=ls())
library(MASS)
library(foreign)
library(Matching)
library(ebal)

### simulation parameters
# treated # control  # no of vars
tr <- 50; co <- 100; nvars <- 3
## Pick Design 
# Design A: Equal unit variances in both groups. 
# Design B: Unequal variances: 1:5 in treatment and :5 in control group. 
# Design C: Equal unit variances. Squared terms of X are included for all estimators, but not outcome
Design <- "C"
# Covariance Matrix
Sigma.tr <- Sigma.co <- diag(rep(1,nvars))
if (Design == "B") {
Sigma.tr <- diag(c(rep(1.5,nvars)))
Sigma.co <- diag(c(rep(.5,nvars)))
}
# mu treated and controls
mu.x.tr <- rep(0,nvars)
mu.x.co <- rep(.2,nvars)
# beta for linear link
betas <- rep(1,nvars)
# number of sims
sims <- 1000
# storage for different estimators
mods  <- c("RAW","PS","MD","PSMD","GM","WPS","EB")
store <- matrix(NA,sims,length(mods))
colnames(store) <- mods

### MC experiment starts here
for(i in 1:sims){

# simulate covars
 # Covars treated
  x.tr <- mvrnorm(tr,mu=mu.x.tr,Sigma=Sigma.tr, empirical = F)
 # Covars controls
  x.co <- mvrnorm(co,mu=mu.x.co,Sigma=Sigma.co, empirical = F)

# treatment indicator
W <- c(rep(1,tr),rep(0,co))

# combine data
X <- rbind(x.tr,x.co)

# simulate outcome
e <- rnorm(nrow(X),0,.5^2)
Y <- X%*%betas + e

# implement estimators
if (Design == "C") {
X <- cbind(X,X^2)
}

# raw
store[i,"RAW"] <- mean(Y[W==1]) - mean(Y[W==0])
# PS dist
PS.out <- glm(W~X,family=binomial(link = "probit"))
PS <- PS.out$linear
store[i,"PS"] <- Match(Y=Y,Tr=W,X=PS,estimand="ATT",M=1,BiasAdjust = FALSE)$est
# MD
store[i,"MD"] <- Match(Y=Y,Tr=W,X=X,estimand="ATT",M=1,BiasAdjust = FALSE)$est
# MD PS combined (with othorgonalized X)
L <- cbind(rep(1,length(PS)),PS)
L <- X-(L%*%(solve(t(L)%*%L)%*%t(L)%*%X))
Xortho <- cbind(PS,L)
store[i,"PSMD"] <- Match(Y=Y,Tr=W,X=Xortho,estimand="ATT",M=1,BiasAdjust = FALSE)$est

# GenMatching
suppressWarnings(gout          <- GenMatch(Tr=W, X=X, BalanceMatrix=X, estimand="ATT", M=1, weights=NULL,
                          print.level=0))
store[i,"GM"] <-  Match(Y,W,X,BiasAdjust=FALSE,Weight.matrix=gout,estimand="ATT",M=1)$est

# weighting on PS
PS.pr <- PS.out$fitted
d <- PS.pr/(1-PS.pr)
Y.co.est <- sum(Y[W==0]*d[W==0])/sum(d[W==0])
store[i,"WPS"] <- mean(Y[W==1])-Y.co.est

# Entropy Balancing
out.eb <- ebalance(
              Treatment=W,
              X=X,
              print.level = -1
              )
d <- out.eb$w
store[i,"EB"] <- mean(Y[W==1]) - (sum(Y[W==0]*d)/sum(d)) 

# output
if((i%%100)==0){
cat("MC Sims to go",sims-i,"\n")
}

}

# Summarize
Bias  <- apply(store,2,mean)
MSE   <- apply((store-0)^2,2,mean)
out <- rbind(Bias*100,MSE*100)
out <- rbind(out,abs(out[1,]/out[1,ncol(out)]))
out <- rbind(out,abs(out[2,]/out[2,ncol(out)]))
rownames(out) <- c("Bias","MSE","Bias / Bias EB","MSE / MSE EB")
round(out,2)
